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An important topic of interest in imaging is the construction of protocols that are not diffraction limited. This can 
be achieved in a variety of ways, including classical superresolution techniques or quantum entanglement-based 
protocols. Here, we consider superresolving imaging in the far field using higher-order intensity correlations. 
We show that third and fourth order correlations can improve upon the first and second order correlations that 
are traditionally used in classical optics and Hanbury Brown and Twiss type experiments. The improvement is 
achieved entirely by post-processing of the data. As a demonstrator, we simulate the far field intensity distribution 
of a circular aperture that emits thermal light and use maximum likelihood estimation to determine the radius 
of the aperture. We compare the achieved precision to the Cramer-Rao lower bound and find that the variance 
of measurements for the third and fourth order correlation functions are indeed closer to the Cramer-Rao bound 
than that of the second order correlation function. The method presented here is general, and can be used for all 
kinds of incoherent emitters, geometries, and types of noise. 

PACS numbers: 42.30.-d, 42.50.St, 03.67.-a, 03.65.Ta 


I. INTRODUCTION 

Visualization of natural phenomena has played a central role 
in the development of science and technology. Indeed, the in¬ 
vention of the telescope shed light on the motion of planets and 
stars and the conception of the microscope allowed investiga¬ 
tions into the ingredients of life at the microscopic level. His¬ 
torically, every time a new imaging technique was introduced, 
science has leaped a step forward. Most recent advances in 
imaging include exoplanet detection [ 1 ] and the velocity mea¬ 
surement of molecular markers along DNA [2]. However, the 
wave nature of light dictates that there are physical limits to 
the resolution of optical telescopes and microscopes, as for¬ 
mulated by Abbe in 1873 [3]. In order to see smaller details 
in microscopy we could illuminate with light of shorter and 
shorter wavelengths, but this is not always practical: highly 
energetic light may destroy biological samples, and in astron¬ 
omy shorter wavelengths are increasingly difficult to access. 
It is thus useful to find alternative techniques to improve the 
resolution of imaging methods that overcome the diffraction 
limit. Imaging techniques that yield finer details than that dic¬ 
tated by the Abbe limit are referred to as superresolving tech¬ 
niques. In microscopy techniques such as photo-activated lo¬ 
calised microscopy (palm) [4], stochastic optical reconstruc¬ 
tion microscopy (storm) [5] and stimulated-emission deple¬ 
tion microscopy (sted) [6] achieve superresolved images using 
fluorescent markers. These methods combine standard inten¬ 
sity measurements with non-linear effects or stochastical pro¬ 
cesses in combination with prior information about the sample 
preparation and post-processing to achieve the pursued goal of 
superresolution. 

Alternatively, superresolving imaging may also be achieved 
using higher-order intensity interferometry [7-9]. Hanbury 
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Brown and Twiss (HBT) in their seminal experiments demon¬ 
strated that the second order intensity correlation function is 
proportional to the Fourier transform of the intensity distri¬ 
bution of a thermal light source [10]. Measuring the intensity- 
intensity correlations therefore provides a method to access the 
spatial distribution of light emitters. The question then arises 
how and to what extent higher-order intensity correlations can 
be used to improve the resolution in imaging. Early propos- 
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FIG. 1. (a) Spatially incoherent, planar source S emitting into 

the far field. The intensity measured at the positions ri,..., r„ 
can be correlated to give the n-point intensity correlation function 
... ,r„), which contains information about the geometry 
of S. (b) Implementation simulated here: a circular aperture of ra¬ 
dius 100 pm emits monochromatic thermal light with a wavelength 
of 633 nm, which is recorded with a CCD camera. The intensity mea¬ 
sured by the different pixels of the CCD can be correlated to calculate 
G(")(ri,...,r„) as in (a). The task is to estimate the radius of the 
aperture from G*"^(ri,... ,r„). 
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als using entangled -photon states promised an increase in 
resolution cx 0{N) [11]. Recent more practical multi-photon 
correlation techniques achieve a similar scaling [7-9], relying 
on the precise estimation of certain parameters. To ensure the 
optimal performance of such techniques it is thus necessary to 
obtain the best estimates possible. In this paper for the first 
time rigorous estimation theory is used to demonstrate how 
higher-order intensity correlation measurements can yield a 
resolution improvement over conventional second order inten¬ 
sity interferometry. The general setup is shown in Fig. 1 . As a 
demonstrator we investigate a planar circular source that emits 
spatially incoherent monochromatic thermal light, which is 
recorded in the far field by n detectors at ri,..., r„. The 
recorded intensities are then correlated in order to produce the 
n‘*'-order intensity correlation function. In this paper, we pro¬ 
vide a general framework that allows us to investigate intensity 
correlations of arbitrary order and discuss the advantages and 
disadvantages with respect to improved resolution produced by 
higher-order intensity correlation measurements. We present 
a rigorous analysis of the HBT set up, explicitly calculating 
the corresponding likelihood function and Fisher information, 
allowing us to determine the lowest possible variance in esti¬ 
mates of the source. We also perform a maximum likelihood 
estimation, enabling us to achieve this bound. An integral part 
of our procedure is to treat unknown quantities as parameters 
to be estimated, every unknown therefore increasing the size of 
the estimation problem. We then demonstrate the method with 
a concrete example, making use of simulated data, to obtain 
precision estimates of the radius of a thermal source. We also 
show how these results can be used to estimate the dimensions 
of any source, regardless of its geometry or its photon statis¬ 
tics, as long as the total photon number is not deterministic (i.e. 
An > 0). For any such source, our method allows us to cal¬ 
culate the Fisher information, providing maximum likelihood 
estimates of the source dimensions via the method of scoring. 
By explicitly calculating the Fisher information for different 
correlation orders n, we then show that in some instances the 
Fisher information increases with correlation order. The exact 
relationship between the Fisher information and the correla¬ 
tion order depends strongly on the particular arrangement of 
the detectors and the source geometry. 

The idea of exploiting spatial intensity correlations of or¬ 
der n > 2 to improve resolution has been suggested previ¬ 
ously [7, 8, 12]. In particular in the context of ghost imag¬ 
ing [13-17], the idea of exploiting higher correlation orders 
has received great attention [18-22]. However, to date there 
has not been a rigorous theoretical approach to quantify the 
possible enhancement in resolution when using higher-order 
correlation measurements for HBT-type experiments. We pro¬ 
vide that explanation here and demonstrate its use. The effects 
of photon losses are included in our model and we show that 
higher-order correlations can continue to perform well, even 
in the presence of loss. 

The paper is organized as follows: In Section II we present 
the theoretical model used to determine the higher-order cor¬ 
relations of the intensity produced in the far field by the inves¬ 
tigated source. We present the exact mathematical expression 
for the corresponding n‘’'-order intensity correlation function 


and show that in the case of a thermal source it can be written 
in terms of the permanent of the correlation matrix involving 
the complex degree of coherence. In Section III we present 
the theory of parameter estimation, which involves the use of 
optimized estimators that take the measured data to return op¬ 
timal estimates of the parameters of interest. This allows us to 
evaluate their performance by use of the Cramer-Rao bound, 
expressed in terms of the Fisher information, i.e., it enables us 
to derive an explicit expression for the lower bound of the vari¬ 
ance in the estimates of the considered source. In Section IV 
we present an explicit protocol for obtaining the data of an n***- 
order intensity correlation measurement, incorporating statis¬ 
tical uncertainties as well as additional noise due to, e.g., the 
limited detection efficiency of the detectors. This enables us to 
calculate the covariance matrix for the given problem, which 
in turn gives access to the Fisher information, the Cramer-Rao 
bound, and the maximum likelihood estimation procedure. In 
Section V we present detailed numerical simulations for dif¬ 
ferent source geometries and calculate the Cramer-Rao bound 
for the source dimension for different orders of the intensity 
correlation function, taking into account two kinds of noise, a 
constant noise factor at each detector and Gaussian distributed 
noise due to detector losses. Finally, in Section VI we discuss 
our results and present our conclusions. 


II. N-POINT INTENSITY CORRELATION EUNCTIONS 

We consider a semi-planar source S emitting monochromatic 
spatially incoherent thermal radiation that is observed in the 
far field. The setup is shown in Fig. la. The second order 
intensity interference for experiments that measure the equal¬ 
time two-point intensity correlation reads 

G(2)(ri,r2) = {■.E^-'>{r,)E^+\r,)E^-\r2)E‘^+\r2):) 

cx (:a^(ri)a(ri)o’''(r 2 )a(r 2 ):) , (1) 

where E^^'> are the positive and negative frequency parts of 
the electric field, a(r) and a^(r) are the annihilation and cre¬ 
ation operators of the field at position r, :: denotes normal or¬ 
dering, and Vi is the position of the i**' detector in the far field. 
When viewed in the far field of the source, interference fringes 
can be observed in r 2 ) under certain conditions [10]. 

However, we can also consider the equal-time n-point inten¬ 
sity correlation G*^"^ cx (clir^i increasing 

the correlation order n can lead to an increased visibility of 
the interference fringes, suggesting we may be able to extract 
more information from the higher correlation orders [8, 23]. 

To obtain the higher-order intensity correlations we consider 
the particular measurement scheme displayed in Fig. 2. An 
array of pixels arranged along an axis parallel to the surface of 
the source allows to measure the intensity at a discrete number 
of positions. The benefit of such a setup is that we can capture a 
large amount of data in one frame. In this section we establish 
the theoretical model of the intensity correlations to any order; 
the estimation procedure will be discussed in detail in Sections 
III and IV. 
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Thermal light exhibits a Gaussian zero mean P-representation. 
We can therefore apply the Gaussian moment theorem [25] and 
make the simplihcation 


crGSri 1 

where S'„ is the symmetric group containing all permutations 
of the set n}. This allows us to write as 

n 

G(”)(ri,... ,r„) = ^ ]J(a’l'(ri)a(r„(i))), (5) 

(y^Sn 1 


where we have defined = Ka. We can make a fur¬ 

ther simplification by introducing the complex degree of co¬ 
herence, defined as [25] [26] 


7(i'i,r2) 


_(a^(ri)a(r2))_ 

[(at(ri)a(ri))(at(r2)a(r2))]^/^ 


( 6 ) 


This allows us to write the n-point intensity correlation func¬ 
tion as 


n 

G(")(ri,...,r„) = ]J(a'l'(ri)a(r,))7(r,,r^(i)), (7) 

o-GS^ i=l 


FIG. 2. Measuring the higher-order intensity correlation functions 
with an array of pixels, (a) The second order intensity correlation at 
pixel Xi — 7 is calculated as the correlation between the intensities 
at a fixed pixel S 2 = 13 (shown as a darker pixel) and at the pixel 
Xi = 7. (b, c) The third order correlation is defined using two detec¬ 
tion schemes, namely via a single fixed pixel that is correlated twice 
(detection scheme 1), or two fixed pixels (detection scheme 2). (d, 
e) The fourth order intensity correlation is defined analogous to the 
third order. The reference pixel separation d is a dimensionless num¬ 
ber, and the center-to-center separation of adjacent pixels is taken as 
5.3 pm throughout the paper. 


To calculate the n-point intensity correlation we make use of 
the optical equivalence theorem [24], which states that the ex¬ 
pectation of a normally ordered product of creation and annihi¬ 
lation operators can be replaced by their left and right eigenval¬ 
ues, respectively, if the expectation is replaced by an ensemble 
average weighted by the P-representation of the state. That is, 
mathematically, we can write 

(/(a^a)) = J P{a)f{a*,a)d'^a = {f{a*,a))p, (2) 

where / is any normally ordered function of the creation and 
annihilation operators, the first expectation is the quantum me¬ 
chanical average, and the subscript P on the second expecta¬ 
tion signifies that it is an ensemble average taken with respect 
to the quasi-probability distribution P. With the help of the 
optical equivalence theorem, the n-point intensity correlation 
G(")(ri,...,r„) can thus be written as 


( 3 ) 


where we omitted the constant of proportionality K for brevity. 
From this expression we see that the n-point correlation for 
Gaussian light is equal to the permanent of a matrix F 


G(")(ri,...,r„) =Perm(r), (8) 

where 

Tij = [(a1'(r,)a(G))(a’^(rj)a(rj))]^^%(ri,rj). (9) 

In general the permanent of a matrix is difficult to calculate. 
Therefore, for larger n it is increasingly costly to calculate the 
n-point correlations, this may be the main limiting factor for 
the use of higher-order correlation functions in imaging. 

The advantage of writing G*-"^ (ri,..., r„) in terms of the 
complex degree of coherence is that the complex degree of 
coherence in the far field paraxial regime is given by the two- 
dimensional Fourier transform of the intensity distribution of 
the source [25] . This result is known as the Van Cittert-Zernike 
theorem [27, 28]. We can therefore calculate the far field in¬ 
tensity correlations of any order for any source geometry, pro¬ 
vided that we can determine the Fourier transform of the in¬ 
tensity distribution. 

As an example, consider a circular source of uniform inten¬ 
sity and angular diameter = 2tan“^(a/L) « 2a/L where o 
is the radius of the source and L is the distance from the source 
to the observation plane. The far field complex degree of co¬ 
herence is the two-dimensional Fourier transform of a circle 
with radius a [29] 




2Ji - r2|) 

(ii9fc|ri - ral) 


(:nr=ia’^(G)a(G):) = (nr=i«*{G)a(i'i))p • 


( 10 ) 
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HG. 3. (Color online) The second (dashed), third (dotted) and fourth 
(solid) order intensity correlation functions of Eq. (11) as a function 
of the separation between the scanning pixel r i = x and the reference 
pixels at S 2 for the remaining arguments, with normalised intensity 
(a^(ri)d(ri)) = 1 and 'd = 5 x 10~* rad. The width of the curves 
is directly proportional to the angular diameter of the source t?. The 
higher visibilities of and G*-^^ over G*'^^ suggest that higher- 
order correlation functions may outperform regular Hanbury Brown 
and Twiss estimation of the source diameter. In this paper, we show 
rigorously using estimation theory that this can indeed be the case 
and specify the conditions. 


where Ji is the first order Bessel function of the first kind and 
k is the wavenumber. The n-point correlation function then 
becomes 




•■,fn) = 


creSr^ i=l 


2Ji {^dk\r, - r^(i)l) 
(I'dklri - r<^(,)|) 


( 11 ) 


Fig. 3 shows the second, third and fourth order intensity cor¬ 
relation functions for a uniform disc along a one-dimensional 
detection device. 


III. ESTIMATION THEORY 


In the previous section we saw that intensity correlation mea¬ 
surements in the far field depend on the parameters describing 
the source geometry. In this section we examine how we can 
use these measurements to practically obtain spatial informa¬ 
tion about the source. To this end we employ parameter esti¬ 
mation theory, which involves the use of an estimator 6 that 
takes the measured data x = (xi,..., Xm) in M measure¬ 
ments (the pixels in the detector) and returns estimates of the 
parameters of interest 0 = (0i,..., 0;). In order to extract the 
spatial information in the most efficient way, we can apply a 
maximum likelihood estimation procedure. Maximum likeli¬ 
hood estimation relies on maximisation of the joint probability 
distribution of our data, and we therefore need to characterise 
the probability distribution from which the correlation func¬ 
tions are sampled [30]. 


Once we have determined the conditional probability 
p(x|0) of obtaining the measurement outcomes x given the 
values of the parameters 6 = (0i,..., 0;), we can determine 
the performance of our estimates using the Cramer-Rao bound 
(CRB). The CRB provides a lower bound for the variance of 
our estimators in terms of the Fisher information matrix 

Var(0,) > [I(0 )]t1, (12) 

where the Fisher information matrix 1(0) is given by 

91n[p(x|0)] \ / ain[p(x|0)] 

36, ) \ 39, 

(13) 


X 


Note that the sum over x may be an integral if the values of 
Xi form a continuum. In practice, intensity measurements in 
modern optical detectors yield a digital signal with discrete 
values. 

The exact probability distribution p(x|0) of the correlation 
functions will be a complicated expression depending on the 
number of images N —not to be confused with the number of 
pixels M in each image. However, since the measurements 
of the correlation functions are averages over a (preferably) 
large data set, the central limit theorem dictates that these 
measurements will be normally distributed [31]. Assuming 
that we make N measurements of the correlation functions 
at M discrete detector positions, the data will follow an M- 
dimensional normal distribution: 

exp(-i(x-p(0))TC-i(0)(x-p(0))) 

p(x 0 = -^^- , - -^ , 

v/(2^)^|C(0)| 

(14) 


where p is the tuple of expectation values of the distribu¬ 
tion at each of the sampling points: p"'" = ((xi),..., (xm)), 
^ denotes the transpose, and C is the covariance matrix be¬ 
tween pairs of measurements = Cov(xi,Xy) = (x^xy) — 

(x*)(Xj). 

For a multivariate-normal distribution the elements of the 
Fisher information matrix are given by [30] 


mh 




+ -Tr 

[c-^c- 

,dC' 

\3dJ 

K36jj 

2 

[ 36, 



[Ii(0)],, + [l2(0)],„ (15) 


where we define the first term (depending on p) as [li(0)]iy, 
and the second term as [l2(0)]iy. Since the optical field ex¬ 
hibits strong transverse correlations in the detection plane, the 
covariances will not be negligible, and we must therefore ex¬ 
plicitly evaluate these covariances in order to perform max¬ 
imum likelihood estimation. In the next section we explain 
how the correlation functions are measured and incorporated 
into the parameter estimation procedure. 

Finally, in order to find the maximum likelihood estimate 
we use an iterative method called scoring [30]. The process is 
described by the recursion relation 


l(0(fe))0(fc-ti) ^^g,(fc))g,(fe) 


c)ln[p(x|0)] 


30 


6)=0('=) 


(16) 
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where 9^^'^ is the iteration of the parameters 8, and 
X\g^g{k) denotes evaluation of the quantity X at the value 
9 = 9^^\ In order to begin the scoring algorithm we require 
an initial value 9^'^\ Provided the initial value is sufficiently 
close to the actual value, the algorithm should continue with¬ 
out difficulty. If no prior knowledge exists about the param¬ 
eters to be estimated, approximate values can be obtained by 
simple methods (which by no means achieve the CRB) that can 
then be used as the initial values 0^°^. 


IV. MEASURING THE CORRELATION EUNCTIONS 

For simplicity we suppress the y dependence in = {xi, yi) 
and consider only the one-dimensional problem where a sin¬ 
gle “moving” detector Xi scans across a set of M discrete 
positions xi,..., xm and the remaining n — 1 detectors are 
kept fixed (see Fig. 2b-e). We refer to the fixed detectors as 
the reference pixels and write the reference pixel positions as 
X 2 = S 2 ,---,Xn = Sn- We distinguish between two detec¬ 
tion schemes, namely one where all reference pixels are iden¬ 
tical (scheme 1), and one where all reference pixels are dif¬ 
ferent (scheme 2). Taking N images means performing N in¬ 
dependent measurements of the intensity I{xi) at all pixels 
Xi = Xi,..., Xm and calculating the sample average of the 
intensity moments 

^ N n 

(17) 

^ k=li=l 

where Ik{xi) is the fc*'' measurement of the intensity at posi¬ 
tion Xi. Given the reference pixels {s 2 ,..., s„}, we can ab¬ 
breviate ©(") (xi, S 2 , ■ ■ ■, Sn) = and the data that 

is used in the estimation procedure is then given by x = 
..., This is still quite a general de¬ 

scription and includes, for example, the experimental arrange¬ 
ment used in Ref. [8]. We use 0^”^ to denote a measurement 
of the correlation function. This is not to be confused with 
the true correlation function as given by Eq. (7). It is an 
important distinction since the measured correlation 0(”)(a;i) 
is a random variable due to the finite size of the sample N, 
whereas G^'^\xi) is the expectation value of the correlation 
function, only in the limit N ^ oo do the two coincide. The 
relation between the standard statistical quantities and the cor¬ 
relation functions are collated in Table I. 

In addition to the statistical noise due to the finite sample 
size N, any realisable detection scheme will introduce addi¬ 
tional noise into the measurements. One important source of 
noise is reduced detection efficiency of the pixels. Often this 
is treated as a constant parameter p. However, when calculat¬ 
ing the intensity correlations we necessarily sample the higher 
moments of the detector noise. It is therefore important that 
we acknowledge the random nature of the noise in order to 
correctly deduce its effects. Physically we would expect the 
noise to be sharply peaked around some constant value with 
some small but non-zero variance. We would also expect the 


random noise to be independent across the pixel array since the 
pixels themselves are independent. We model this additional 
noise as uncorrelated, normally distributed noise with mean 
and variance {r]{xi)) = v and {ri{xi)‘^) — {rj^Xi))^ = re¬ 
spectively [32]. We can therefore write 

^k{Xi^ — rjk{Xi]ik{Xi^ t (18) 

where Ik{xi) is the fc**' realisation of the random intensity at 
pixel Xi as measured by an ideal detector, and rjkixi) is the fc**' 
realisation of the noise at pixel Xi. The expectation of 0^”! {xi) 
is then given by 

fii = {<3^'^\x^,S2 ■ . ■ ,Sn)) (19) 

1 ^ 

= '^{Ikixi)Ikis 2 ) . . . IkiSn)) 

^ k^l 

I ^ . 

= '^{ik{Xi)ik[s 2 ) ■ ■ . 4(s„)) 

k^l 

X {r]kixi)r]kis2) ■ ..rjkisn)) , 

which is to be used in Eq. (14). We assume that the noise 
and intensity are stationary random variables, so we can im¬ 
mediately perform the sum removing the factor N~^. The 
first factor in Eq. (19) is, by definition, the intensity correla¬ 
tion {ik(.Xi)ik{s2) ■ ■ ■iki.Sn)) = G(”1(xi,S2,---,s„). The 
second factor in Eq. (19) is in general some combination of 
moments of the normal distribution characterised by v and ^ 
(see Appendix). 

As shown in Eig. 2a, we calculate the two-point intensity 
correlation 0^^! (a;, S 2 ) between any pair of pixels (a;, S 2 ) as a 
function of the position of the pixel at x and the stationary posi¬ 
tion of the second pixel at 32- Note, however, that when mea¬ 
suring a correlation function in this way the individual data 
points that are calculated may not be independent; The corre¬ 
lation between any pair of data points, e.g., 0l^)(a:i, S 2 ) and 
0^^1 {x 2 ^ S 2 ), depends on the correlations between all the mea¬ 
sured intensities Ik{xi), Ik{x2) and Ik{s2)\ since the mea¬ 
surement relies on the statistical dependence of these intensi¬ 
ties, the resulting values of iS^^\xi, S2) and 0*^^^ (x2, S2) will 
not, in general, be statistically independent. The same argu¬ 
ment holds for higher-order correlations. In principle it is pos¬ 
sible to avoid correlations between the data points altogether 
by taking each measurement in a completely independent man¬ 
ner. Hanbury Brown and Twiss did just this in their original 
experiments. However, the price paid for taking data in such 
a way is a much greater total measurement time. A more effi¬ 
cient way to collect the data would be to make use of an array 
of detectors (usually pixels of a CCD camera) to take all the 
measurements simultaneously. As long as we are careful to 
take into account the correlations that arise in the data when 
measured in this manner then we are free to use this efficient 
method of data collection. In the following we assume that the 
data is collected in such a way and are careful to calculate the 
correlations in the data explicitly. 

The elements of the covariance matrix in Eq. (14) are given 
by 
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statistics 


n*’'-order intensity correlation measurements 


random variable 
fc*** measurement 
sample mean 
first moment 
covariance 


j^(k) 

m = (X.) 

Cov{Xi,X,) = j^Cov{Xi,Xj) 


I{Xi)I{s2) . . . I{Sn) 

Ik{xi) 7fc(s2) • • • h{s„) 

&^'^\xi) = ^ Ik{Xi)Ik{s2) ■ ■ . 4(Sn) 

(0'") (*0) = (Hxi) I{S2) . . . /(S„)) OC G'") iXi,S2,...,S„) 

Cov (^&^"'>{xi), (S^”-\xj)'^ = ^ Cov {I(xi) I{S 2 ) ■ ■ ■ I{s„),I{Xj) I{s 2 ) ■ ■ ■ 


TABLE I. Statistical quantities and their counterparts in the n'^’'-order intensity correlation measurements. The random variable Xi is the 
product of the intensity measurements at positions Xi, S 2 ,-■ ■, Sn, and the index i £ {1,..., M} runs over all the detector positions (pixels). 
After N measurements, we define a sample mean Xi that is itself a fluctuating quantity. This is not to be confused with the first moment 
(Xi), which is not a random variable. Since the sample mean is an unbiased estimate of the first moment, the expectation value of the sample 
mean is equal to the first moment. The covariance of any pair of sample averages is not equal to the covariance of the variables but N~^ times 
the covariance. This quantifies the fact that taking more data (increasing N) reduces the variation of the sample averages. As N ^ oo the 
sample averages coincide with the expectation values and the Xi are no longer random variables. The first moment is only proportional to the 
correlation function due to the efficiency factor in the measured intensity in Eq. (18). 


■'ij 


1 ^ 

jy-2 ^ ^ {Xi)Ik ('^2) ■ ■ ■ {Xj )7;(s2) ■ . • 5 


( 20 ) 


where Cov denotes the covariance between the measured in¬ 
tensity correlations at pixel Xi and Xj. Since k and I label 


Cjy — ^2 I 'y ^^ {kk(,Xi^Ik{s2^ ■ • • ^k{^^n}^k{Xj^Ik(^^2) • ■ • ^k^^ 

\k=l 
- fJ'ifJ-j. 


Since we treat the intensities as stationary random variables, 
we can simply perform the sums over k and I to obtain 

Cij = ^ (^{I{Xi)I{Xj)I{s 2 f ■ ..I{Snf) - (22) 

— ^ (^Xif Xj , S2j S 2 , . . . 5 Sjif Sji^ 

X {r]{x,)r]{Xj)r]{s2)v{s2) ■ ■ ■ v{sn)v{sn)) - k'ik-j] ■ 

We see that the covariances between our data and 

®(")(xy) depend on correlation functions of order 2n. The 
term {r]{xi)r]{xj)T]{s2)r]{s2) ■ ■ ■ r]{sn)r]{sn))) is evaluated in 
the Appendix. We now have a complete characterisation of 
the probability distribution p(x|0) and can therefore calculate 
the Fisher information to determine the lower bound on the 
variance of our estimates of 6 via the CRB and can also per¬ 
form a maximum likelihood estimation procedure to estimate 


the individual images that are statistically independent, we can 
split the sum into two parts 


)) + ^ {h{Xi)Ik{s2) . . ■Ik{Sn)){h{Xj)Il{s2) ■ . . Il{Sn)) | 
k^l 

( 21 ) 


the dimensions of the source. In the next section we present 
numerical simulations of this estimation procedure. 


V. NUMERICAL SIMULATIONS 

To determine the performance of our estimates, we produce 
simulations of the experiment shown in Fig. lb, where a cir¬ 
cular aperture of radius of a = 100 pm emits uniform thermal 
light of wavelength A = 633 nm. The complex degree of co¬ 
herence of such a source in the far field is given by Eq. (10). 
To simulate the experiment we again make use of the opti¬ 
cal equivalence theorem and the P-representation. The P- 
representation for thermal light takes the form of a complex 
multivariate normal distribution. The simulation of the exper¬ 
iment is then performed by sampling from a 2M-dimensional 
normal distribution corresponding to the real and imaginary 
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parts of a{xi ),..., a{xM), i-S-, Oi{x) = a{x) + ib{x) with 
covariances 

{a{x,)a{xj)) = {b{xi)b{xj)) = ^(7) -f{xi,Xj) 
{a{xi)b{xj)) = 0 . (23) 

Here we have assumed a uniform far field intensity distribution 
of the thermal source, i.e., {I{xi)) = {I{xj)) = (7). 

We include the effect of pixel noise by adding an additional 
normal random variable to each of the intensities with mean 
1 / and variance The intensity correlations are then calcu¬ 
lated from the simulated field by means of Eq. (17), which are 
in turn used to estimate the dimensions of the source param¬ 
eters. We require averaging over a large set of data in order 
to apply the central limit theorem and treat the data as nor¬ 
mally distributed. Here the parameter (7) is also treated as an 
unknown parameter to be estimated. The importance of this 
cannot be overstated; if we instead treat (7) as a constant, any 
slight deviation from the exact value can lead to catastrophic 
failure of the estimation procedure. It is therefore imperative 
that the unknown parameter (7) should be treated as a nui¬ 
sance parameter in order to perform the maximum likelihood 
estimation. Although there appear to be four unknown param¬ 
eters: = (a, (7), ly, <;), we will see in Section V A and V B that 
we can often combine the parameters (7) and i/ into the new 
parameter (les) = 1^(7), ly and <; can be combined into the 
new parameter x = These represent the effective inten¬ 
sity recorded by the detector in the presence of inefficiencies 
characterised by (rj) = ty and the ratio of the average pixel 
inefficiencies to the standard deviation of the inefficiencies. 

Once we generated the simulated data from the 2M- 
dimensional normal distribution described above and added 
the noise, we estimated the parameters a, (les) and x based 
on different orders of intensity correlations and the scoring 
method. The maximum likelihood estimation procedure was 
repeated 1000 times such that a statistical variance (Aa)g;jjj 
for the simulated data could be calculated, and we can com¬ 
pare this to the lower bound on the variance (Aa)Qj^B based 
on the Cramer-Rao bound; 

(^o)cRB ^ [I]oo^ : (24) 

where we now must find the inverse of the 3x3 Fisher infor¬ 
mation matrix I. There are two main cases to consider. 


A. Constant detector loss 

First, we analyse the special case of a detection system 
with constant loss for each pixel (i.e., = 0). For this 


particular noise model, the choice of reference pixel posi¬ 
tions does not affect the noise terms in Eqs. (19-22), since 
{Vk{xi)T]k{s 2 ) ■ ■ ■ rjkisn)) — for all choices of reference 
pixel positions. We find the choice S 2 = S 3 = • • • = s„ = 
s = [M/ 2 J to be of particular interest since it simply involves 
taking powers of the measured intensity 7(s)"“^, and it is the 
central pixel on a one-dimensional CCD. This allows us to 
compare the effects of the post processing without the need to 

2 3 4 2 3 4 5 



2 3 4 2 3 4 5 


Correlation order Correlation order 


FIG. 4. (Color online) Results of the numerical estimation of the 
source diameter a (pm) and the corresponding variance (Aa)^, in 
comparison with the Cramer-Rao bound (CRB). Choosing the refer¬ 
ence pixels in a dishibuted manner as in Fig 2c, e leads to a CRB 
value of the variance as shown in (a), with the estimate and standard 
deviation shown in (b). Here we have chosen a detection efficiency 
V = Q.b with ? = 0.01. Due to the computational complexity of the 
problem only correlation orders up to n = 4 have been calculated. 
Choosing the central pixel as the reference pixel as in Fig 2b, d leads 
to a CRB value of the variance as shown in (c), with the estimate and 
standard deviation shown in (d). This configuration does not allow us 
to include the detector efficiency as a random variable in the estima¬ 
tion procedure (see text for details). For comparison, the numerical 
values in this figure are collated in Tables II and III. 


consider complications regarding the exact placement of the 
reference pixels, S 2 , ■ ■ ■ ,Sn- The inherent simplicity of this 
arrangement also allows us to calculate the correlations up to 
arbitrary order since the correlation functions and G^^"^ 
take the compact forms 


G'''^"‘\xi,Xj,s,. 


..,s) 

..,s) 


= (7)"(n-l)![l + (n-l)|7(x.,s)|2], (25) 

= 2)!{1 + \x{x^,Xj)\‘^ 

+ {n-2) [ 2 Re( 7 (a;,,a;y) 7 (s,a;i) 7 (a:y,s)) -F | 7 (a:i,s )|2 -F \x{xj,s)\^ + 17 ( 0 ;^, s)H 7 (xy, s)p] }. 
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FIG. 5. (Color online): Variance (Aa)^ of the estimator a for 
against the standard deviation of the noise (dimensionless). Average 
detector efficiencies v = 0.2 (dashed), v = 0.5 (dotted), and v — 0.9 
(solid). 


We can therefore make the re-parameterisation suggested 
above, leaving us with the parameters 6 — {a, (les)) to esti¬ 
mate. Fig. 4c shows the Cramer-Rao lower bound on the vari¬ 
ance for our estimate of a for the four correlation functions 
to and Fig. 4d shows the estimate of a with the ac¬ 
tual standard deviation Aosim- The numerical results are also 
collated in Table II. 

The best estimates of the spatial dimensions of the source 
occur for n = 3, whereas the estimates get progressively worse 
as the correlation order is increased beyond third order. There¬ 
fore, to extract the maximum amount of spatial information 
from our data, correlations of third order should be used. We 
stress the importance of this finding as it requires no additional 
measurements to be made other than those made to measure 
G^^\ Indeed, in principle it would even be possible to use the 
data collected by Hanbury Brown and Twiss to measure the 
angular diameter of Sirius [10] with improved resolution. We 
also note that there is nothing in our treatment that uniquely 
picks out the spatial correlation functions, in the same manner 
we could equally discuss temporal correlations (see, e.g., [9]). 
Interestingly, estimates of the effective intensity (leff) do not 
follow the same pattern as those for a. If we wish to estimate 
(Jeff) the best performance is given by G*-^^, with higher-orders 
performing worse. 


TABLE II. Results of the maximum likelihood estimation for the cor¬ 
relation functions G^^^ to G^®^ and the Cramer-Rao lower bounds. 
The estimation procedure was performed on 1000 simulated data sets. 


n 

a (pm) 

(Aa)iin, (pm^) 

(Aa)cRB (hui^) 

2 

99.978 

0.181 

0.162 

3 

99.968 

0.151 

0.150 

4 

99.932 

0.249 

0.244 

5 

99.826 

0.543 

0.434 



FIG. 6. (Color online): Standard deviation of the estimator a for G^^^ 
(dashed) , G^®^ (dotted) and (solid) as a function of reference 
pixel separation d — |si — Si+i| for detection scheme 2, given a 
circular aperture of radius a = 100 pm. The grey vertical lines cor¬ 
respond to the first and second zeroes of the Bessel function Ji. 

B. Detector loss as a random variable 

Second, we demonstrate the effect of a small non-zero rep¬ 
resenting a system with uncertainty in the detector loss mech¬ 
anism. The effect of this additional noise is shown in Fig. 5, 
where we plot the variance of the estimator for the second order 
intensity correlation function against the standard deviation <j. 
As expected, the addition of noise in the detection process re¬ 
duces the precision in our estimator. 

In order to perform the estimation, first we must evaluate the 
second term in Eq. (19), which is the n**' moment of the noise 
distribution. Having considered the case where all reference 
pixels are the same in the previous section, we now restrict 
ourselves to only considering cases where no two reference 
pixel positions are the same, i.e. S 2 ^ S 3 s„, such 

that we can determine the effect of separating the reference 
pixels. In this regime the n**' moment of the noise distribution, 
{r]{xi)'q{s 2 ) ■. ■ v{sn)), is given by 

n 

{rikixi)r]kis 2 )... rjkisn)) = ^ 

l + ■ ( 26 ) 

i=2 ) 

We therefore find it necessary to re-parameterise the problem 
using the parameters 6 = (a, (Jeff),x) as mentioned above. 
Eq. (26) represents an n**' moment of the noise distribution, 
and the Kronecker deltas arise from the independence of the 
distribution for individual pixels. The 2n*‘^ moment of the 
noise distribution appearing in Eq. (22) is calculated in the Ap¬ 
pendix. 

In order to find the optimum position of the reference pixels, 
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HG. 7. (Color online): Standard deviation of the estimator a for 
(dashed) , G^®^ (dotted) and (solid) as a function of reference 
pixel separation d = |si — Si+i| for detection scheme 2, given a slit 
of width a = 200 pm. The grey vertical lines correspond to the first 
and second zeroes of the sine function. 


we define the dimensionless number d = |si — Si+i|, the sep¬ 
aration between adjacent reference pixels, and plot the stan¬ 
dard deviation as a function of d. Fig. 6 shows the standard 
deviation for to as a function of d. Interestingly, 
the higher-order correlations outperform G*^^^ only for some 
values of d. For G^^^ we find that the optimum positions cor¬ 
respond to separations where the two reference pixels become 
uncorrelated. This occurs whenever the complex degree of co¬ 
herence between the two pixels is approximately equal to zero. 
Since the complex degree of coherence for the system is pro¬ 
portional to the Bessel function Ji, the optimum separations 
d correspond to the zeroes of this function. For G*^^^ the ex¬ 
act position of the optimum is more complicated, due to the 
fact that the zeroes of Ji are not uniformly distributed. How¬ 
ever, the optimum positions are approximately located at the 
position where adjacent reference pixels are uncorrelated from 
their nearest neighbours. Tab. Ill shows the variance of the es¬ 
timators for the first three correlation functions as calculated 
from the CRB and directly measured in the simulations. We 
see that the measured variance in our estimators closely fol¬ 
lows that obtained from the CRB. The results of the maximum 


TABLE III. Results of the maximum likelihood estimation for the 
correlation functions G^^Go G^^^ and the Cramer-Rao lower bounds, 
with a reference pixel separation of d = 182 that corresponds to the 
first zero of Ji. The estimation procedure was performed on 1000 
simulated data sets. 


n 

a (pm) 

(^«)sim (pm^) 

(Ao)crb 

2 

99.976 

0.194 

0.175 

3 

99.980 

0.126 

0.123 

4 

99.957 

0.169 

0.157 


likelihood estimation and the Cramer-Rao bound are given in 
Table III. 

Another interesting feature of Fig. 6 is the ability for G^) 
to outperform G^), a feature that does not occur for fixed ref¬ 
erence pixels. This behaviour is reminiscent of the “magic an¬ 
gles” in Refs. [7, 8], where the detectors had to be placed at 
different specific positions (the magic angles) in order to ob¬ 
tain the {n — l)-fold increased sinusoidal modulation in the 
scanning detector. 

The exact relation between the variance of the estimators 
and the correlation order also depends on the geometry of 
the source. Fig. 7 shows the dependence of the variance as 
a function of the reference pixel separation for a slit of width 
a = 200 |xm. The complex degree of coherence for such a ge¬ 
ometry is given by the sine function. Since the zeros of the sine 
function are uniformly distributed, it is possible to achieve in¬ 
dependence for all the reference pixels simultaneously. Fig. 7 
shows that for the optimal choice of d the estimator for G^) 
outperforms G^) and is about as good as 


VI. DISCUSSION & CONCLUSIONS 

We have discussed the exact role of higher-order intensity cor¬ 
relations with respect to parameter estimation of the intensity 
distribution of thermal sources, and demonstrated that it is 
beneficial to post-process the data in such a way as to measure 
intensity correlations of order n > 2. We have also shown how 
the post processing can be optimised with respect to the place¬ 
ment of the reference pixels, in order to find the most informa¬ 
tive measurements. A major benefit of this method is that it 
does not require particularly elaborate experimental arrange¬ 
ments. Indeed, in certain circumstances it is even possible to 
increase the precision simply by taking powers of the mea¬ 
sured intensities. Since we explicitly account for correlations 
between the data points [see Eqs. (20-22)], all the measure¬ 
ments can be made simultaneously, thus reducing consider¬ 
ably the measurement time required to obtain the data. While 
we have framed the discussion in the context of detector pix¬ 
els on a CCD camera, the same methods apply to any array of 
field detectors, including telescopes. By fully determining the 
probability distribution function (PDF) for measurements of 
intensity correlation functions, including the covariance ma¬ 
trix of the correlated data, we are able to determine the Fisher 
information for such experiments. This allows us to calculate 
the maximum achievable precision via the Cramer-Rao bound 
and also to saturate that bound by performing a maximum like¬ 
lihood estimation. 

The techniques presented here can in principle be used to 
estimate the source dimensions of any object that emits in¬ 
coherent light, including quantum emitters. In many cases 
the Gaussian moment theorem does not apply and extra care 
must be taken in the calculation of the correlation functions 
GG)(a;i , S 2 , ■ ■ ■, Sn)- As long as the measured intensity is a 
random variable then the PDF for the data can be considered 
again as a multivariate normal when averaged over many mea¬ 
surements. Also, the estimated parameters need not be spatial 
parameters of the source: as long as the PDF depends on the 
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parameters to be estimated in a deterministic way our proce¬ 
dure can be used to perform the estimation. The inclusion of a 
realistic noise model allows the effect of detection noise to be 
accounted for, thus allowing the estimation to proceed in the 
presence of noise. 

The experimental implementation of this procedure is not 
exceedingly challenging, as it requires the same kind of tech¬ 
niques as used in the original HBT experiment some 60 years 
ago. A number of critical points must be met though. First, 
when measuring the intensity correlation functions, the inte¬ 
gration time of the detectors must be well below the coherence 
time of the radiation to ensure that each measurement captures 
a single longitudinal mode of the radiation. In addition, the 
area of each detector (pixel) must be much smaller than the 
coherence area of the source such that a single pixel can be 
considered as measuring a single transversal mode. This en¬ 
sures that every detection event samples no more than a sin¬ 
gle mode volume (i.e., the speckle size). Second, our use of 
the central limit theorem implies that we have to use a large 
data set that reduces statistical uncertainties to a minimum and 
provides a good signal to noise ratio. This is a well-known re¬ 
quirement for any maximum likelihood estimation procedure. 
Third, we have assumed uniform mean intensity across the en¬ 
tire CCD. In practice this can be challenging as a natural source 
might have a non-uniform intensity profile, or there might be 
additional spurious interference effects at the detector (e.g., an 
etalon effect due to the protective glass of a CCD camera). Fi¬ 
nally, we have assumed that the detector noise is uncorrelated, 
which means that the pixel efficiencies can be considered ran¬ 
dom and there is no cross-talk between the pixels. The latter 
two requirements can in principle be included in the modelling 
of the experiment, but this comes at the cost of a significantly 
increased complexity of the correlation functions. 

The ability of higher-order intensity correlations to dis¬ 
play more information than lower order correlations is often 
a source of confusion. In fact, the original Hanbury Brown 
and Twiss experiment caused great controversy, and its even¬ 
tual resolution heralded the beginning of quantum optics as 
a mature discipline [33]. The improved estimation capability 
is most clearly demonstrated when the reference pixels are all 
the same S 2 = S 3 = • • • = s„. We can then take the output of 
two photodetectors and simply by taking powers of one of the 
outputs we can gain a more precise estimate of the angular di¬ 
ameter of the source. We could understand this increase in pre¬ 
cision by comparing the measurements of the first and second 
order correlation function as in the original HBT experiment. 
Here we use the same set of data when measuring the first or¬ 
der correlation function, i.e., the intensity distribution in the 
far field, as we use to measure the second order intensity cor¬ 
relation, and yet a measurement of the intensity reveals almost 
no information about the source since the intensity of a ther¬ 
mal source in the far field is constant across x. In contrast, the 
second order intensity correlation function is highly dependent 
on X, which allows for a lensless measurement and therefore 
a precise estimate of the angular diameter of the source (see 
Fig. 3). When considered as another method of post process¬ 
ing the data, it is no more surprising that higher-order correla¬ 
tions outperform the second order correlation than the second 


order correlation outperforming first-order intensity measure¬ 
ments. 

Our method is rather general, and there is nothing in our 
treatment that uniquely picks out the spatial correlation func¬ 
tions, in the same manner we could equally discuss temporal 
correlations. Future work will focus on (multi-mode) squeezed 
light and single photon sources. Since all correlation functions 
can be determined from the same data set, it would be advan¬ 
tageous to combine all of the estimates achieved via different 
n into a single estimate. In order to do this properly we would 
need to know exactly how all of the individual estimates are 
correlated to determine the appropriate weights for the com¬ 
bined estimate and its error. The difficulty in determining the 
correlation is in knowing how the maximum likelihood esti¬ 
mation procedure affects the correlation, if at all. Once this 
is known, it should be possible to determine the weights and 
obtain the final estimate. 

In conclusion, we found that higher-order correlation func¬ 
tions can substantially improve estimation of the parameters 
that characterise the geometry of an incoherent light source. 
As long as the source and the detection system are properly 
modelled, the procedure can be implemented with current 
technology. 
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Appendix A: Moments of the noise distribution 


In this Appendix we evaluate the moments of the noise dis¬ 
tribution for the general case and for the case of a Gaussian 
noise distribution as considered in Section IV. First, we need 
to evaluate the term {r]{xi)r]{s 2 ) ■ ■ ■r]{sn)) which appears in 
Eq. (19). In Section V we calculated this term for the Gaus¬ 
sian noise distribution, here we calculate it for the general case. 
Denoting Ji = ( 77 ( 52 ))... {r]{sn)) and S = {s 2 ,..., s„} we 
find 


{f]{Xi)7j{s2) ■ . .77(S„)) 


{v{xi))Ji 

j 


ii Xi ^ S 
if Xi € S' ’ 


(Al) 


if none of the reference pixels are equal, i.e., S 2 7 ^ S 3 7 ^ • • • 7 ^ 
s„. If instead we use detection scheme 1, where all the refer¬ 
ence pixels are the same, we find 


{v{xi)v{s2r 


{r]{Xt))J2 if Xi 7 ^ 52 

( 77 ( 0 ;^)”) if Xi = 52 ’ 


(A2) 


where J 2 = ( 77 ( 52 )"“^). 

Now we evaluate the second term in Eq. (22), 
{r]{xi)r]{xj)r]{s 2 )'^ ■ ■ ■T]{sn)'^)- Since the noise is treated 
as uncorrelated between the pixels, the expectation value 
factorises into (? 7 (xi)) (? 7 (xj)) ( 77 ( 52 )^) ... ( 77 ( 5 ^)^) if x^ 
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and Xj are not equal to each other or any of the positions S 2 ,..., Sn- However, more generally, the expression is more 

complicated. Denoting J 3 = ( 77 ( 52 )^)... (p(s„)^) we obtain 




{r){xi)ri{xj)T](s2f ■. .riisnf) = < 


J3 = a 


{'n{xiY)Jz = b 

{vjxif) iviXjf) T ^ 

{vixi)){v{xj))J3 = d 




and for detection scheme 1 


{ri{xi)Tj{xj)v{s2) ” ) = < 


{ri{xiY)Ji 

{'n{xi)){r]{xY)Ji 

{ri[xi)){r]{xjY'^-Y 

y{v{x3)){v{xif'^~Y 


if Xi = Xj G S 

if Xi = Xj ^ S 

if Xi Y Xj ; Xi,Xj G S 

if Xi Y Xj ; Xi^Xj ^ S 

if Xi Y Xj', Xi G S; Xj ^ S 

if Xi Y Xj; Xi ^ S; Xj G S 


(A3) 


if Xi = Xj = S2 

if Xi = Xj Y S2 

if Xi Y Xj] x^ Y S2; Xj Y S2 , 

if Xi Y Xj = S2 

if Xj Y Xi = S2 


(A4) 


where J 4 = ( 77 ( 52 )^””^)- For a general noise distribution, 
each term in these piecewise functions can be associated with 
a parameter to be estimated. The use of a general noise distri¬ 
bution can therefore be incorporated in the theory but comes 
at the expense of a greater number of estimation parameters. 
The benefit of using a Gaussian noise model is that there are 
only two additional noise parameters corresponding to the first 
and second moments of the Gaussian distribution, or a combi¬ 
nation of them as in section V, {v,x). 


To give a more visual presentation of the noise correlations, 
we can represent the resulting piecewise function, Eq. (A.3), 


in matrix form M^- = {r]{xi)r]{xj)r]{s 2 Y ■ ■ ■ vYnY) 


M = 


r 

d 

d 

e 

d 

d 

d 

e 

...\ 

d 

b 

d 

e 

d 

d 

d 

e 


d 

d 

b 

e 

d 

d 

d 

e 


e 

e 

e 

a 

e 

e 

e 

c 


d 

d 

d 

e 

b 

d 

d 

e 


d 

d 

d 

e 

d 

b 

d 

e 


d 

d 

d 

e 

d 

d 

b 

e 


e 

e 

e 

c 

e 

e 

e 

a 


L • 








■■ 


(A5) 


where we have made use of the fact that / = e since the noise 
is assumed to be the same for all pixels. 
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